#!/usr/bin/python
from visual import *
from particle import *

# vpython configuration
display(title='Fusion Sim', width=600, height=500, center=(0,0,0), background=color.white, autoscale=False,scale=(0.02,0.02,0.02))


# constants
Q = 1.602e-19 # q > 0 by convention
K = 1/(4*pi*8.854187817*(10**-12))
E_MASS = 9.10938e-31
P_MASS = 1.6726e-27
dt = 0.05
fr = 1/dt



electron = particle(9.10938e-31,-Q,vector(2,0,0),vector(0,.1,0))
#electron = sphere(pos=(2,0,0), color=color.blue)
electron.v = vector(0,0,0)


proton = sphere(pos=(-3,0,0), color=color.yellow)
proton.v = vector(0,0,0)



while 1:
  #rate(fr)

  # calculate forces on particles
  
  # electrostatic force
  R = mag(electron.pos - proton.pos)
  F_coul_pe = K*(Q*(-1*Q))/(R*R)
  F_coul_ep = K*(Q*(-1*Q))/(R*R)  
  
  # gravitational force (just for fun)
  
  F_net_e = F_coul_ep * norm(electron.pos - proton.pos)
  F_net_p = -1*F_coul_pe * norm(electron.pos - proton.pos)
  
  
  # basic newtonian motion, ie: integrate a(t) twice  
  electron.pos = electron.pos + (electron.v*dt) + ((F_net_e/(2*E_MASS))*(dt*dt))
  proton.pos = proton.pos + (proton.v*dt) + ((F_net_p/(2*P_MASS))*(dt*dt))
  
 
  # set new velocity for next dt interval
  electron.v = (electron.v*dt) + ((F_net_e/(2*E_MASS))*(dt*dt))
  proton.v = (proton.v*dt) + ((F_net_e/(2*P_MASS))*(dt*dt))

